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Abstract 

Functional Near- Infrared spectroscopy (fNIRS) is an emerging non- invasive brain computer in- 
terface (BCI) modality that measures changes in haemoglobin concentrations in the cortical 
tissue. To date most fNIRS studies have used standard multiple subject /session dependent 
classifiers for neural signal decoding. Such an approach is preferable because of the available 
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large degree of inter-subject and inter-session variabilities in the acquired data. Thus far, no 
quantitative research has been conducted on these variability issues. In this study, we present 
our methodology to overcome the problem of inter-subject /session variabilities by developing 
classifier functions that maintains good BCI performance regardless of variability in the data. 
We analyzed the fNIRS signals acquired from seven healthy subjects, across eight sessions while 
performing two different cognitive tasks. The proposed classifiers are based on support vector 
machines (with Gaussian kernels) and their extensions to multiple subject /session independent 
feature spaces. We performed a number of extensive experiments, which included developing 
a subject and session dependent classifiers as well as further data independent classifiers. Ex- 
perimental results shows that through the proposed method one can improve the overall BCI 
decoding accuracy. 

Keywords: Brain Computer Interfaces, Functional Near-Infrared Spectroscopy, Inter-subject 
variability, Support Vector Machines, RKHS. 

1. Introduction 

Recent studies have manifested the potential of functional near-infrared spectroscopy (fNIRS) 
as an alternative BCI modality to provide some attractive characteristics over the other stan- 
dard BCI systems [1,2]. The method of fNIRS presents neural information which is related to 
cortical hemodynamics and oxygenation status during functional brain activities through con- 
centration levels of oxy- and deoxy-hemoglobin chromophores of the brain (hereafter denoted as 
oxy-Hb and deoxy-Hb, respectively). The oxy-Hb and deoxy-Hb concentration parameters 
provide important information on subjects' cognitive state and can be mapped into external 
BCI commands [3]. 

To date, most fNIRS based BCI studies have been performed offline in order to decode 
mental task relevant hemodynamic signal responses [4,5, 6-8]. For instance, Sitaram et al. 
demonstrated the decoding of brain hemodynamics arising from right- and left-hand motor 
imagery from five subjects [4]. In their study, the authors achieved nearly 90% classification 
accuracy for motor imagery tasks. In another study, Sassaroli et al. showed successful classifica- 
tion of brain hemodynamics from five mental tasks using a simple k-means algorithm [5] . Other 
offline studies, for instance, have proposed a neural network approach for the classification of 
fNIRS signals [6] of other sophisticated machine learning methods to assess the classifiability 
of fNIRS signal features for a BCI [7]. Most of these studies have produced results in the range 
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of 70%-90% classification accuracy. These initial off-line studies were necessary to establish 
fNIRS as a BCI modality, however, for any practical BCI, it is still imperative to be able to 
operate on-line with the classifier outputs being produced in real-time. In this regard, only 
a few researchers have attempted to use fNIRS signals in on-line settings using different BCI 
experimental paradigms [10-13]. To mention some of the studies, Coyle et al. proposed an 
on-line classifier that compared mean concentration of hemodynamic changes during a motor 
imagery task to a baseline measurement [ 10 ]. Utsugi et al. presented a real-time BCI con- 
trol of a toy train using an online classifier [11]. Their classifiers were trained in a calibration 
mode to generate control commands to activate the train. Another online study was presented 
by Abdelnour and Huppert who used a direct topographic map of oxy-Hb and deoxy-Hb over 
the motor cortex to classify between left and right finger tapping [12]. In their study, online 
classification rates ranged from 68.8% to 93.8% for three subjects, however, finger tapping is a 
physical as opposed to a mental task and studies have shown that imagined finger tapping is 
less discriminable than real finger tapping. Notably, Ayaz et al. used a bar-size control task 
with an update frequency of 2 Hz [13]. Five participants were given up to 120s to increase 
the bar-size to 90% of its maximum size beginning at 50%, but an uncharacteristically long 
task window could mean that the increase in bar-size was due to high-amplitude low-frequency 
oscillations. Most of the mentioned studies exhibit potential in the emerging field of fNIRS 
based BCI research ; nevertheless, there exist some research challenges that need special atten- 
tion. In particular, none of these studies had a chance to address the popular BCI problem, so 
called inter-subject /inter-session variability. This is a major challenge especially in single- 
trial BCI paradigms, whether in fNIRS based BCI or the most advanced BCI systems using 
Electroencephalogram signals [14]. Further, this problem occurs either with invasive or non- 
invasive BCI systems. The inter-session variability that arises in neural signal measurements 
obtained during different recording sessions varies dramatically even when the same subject 
is used. Further, inter-subject variability in signal measurements leads to the development of 
subject-dependent algorithms which are also subject to inter-session variability. The presence 
of this problem strongly degrades the classification accuracy of the BCI system and leads to the 
development of subject /session dependent classifiers [15]. 

In EEG based BCI research, there have been numerous studies focused on overcoming BCI 
variability issues. For instance, there is the work of Matthias Krauledat et al., who first pro- 
posed a zero training method in [16] by extending the applicability of the popular common 
spatial patterns (CSP) algorithm to generalize from one session to other sessions. Their method 
attempts to find prototypical spatial filters from past subject specific sessions, and uses a small 
number of trials of the current session to update the BCI classifier. In another study [17], 
S. Fazli et al. studied an ensemble method which selects a sparse set of subject dependent 
spatio-spectral filters derived from a large database with the recordings of 45 subjects. The 
work showed that a subject-independent classifiers performed almost as well (68% correct) as 
the average subject dependent CSP classifier (70% correct). However, the subject independent 
classifiers' predictions were post-processed with a non-causal bias-correction, which unfortu- 
nately prevents online application. Further examples include the work of F.Lotte et al. who 
in [18] combined different feature extraction and classification methods to compare on their 
ability to discriminate between classes of imaginary movement in unseen subjects. Of all tested 
combinations, a filter-bank CSP classifier that used frequency filters with different bandwidths 
had the best subject independent performance (71% correct). This is actually slightly above 
the subject independent performance of naive log band-power features (68%), and far below 
the best subject dependent classifiers (82%). These three studies indicate that constructing an 
subject independent BCI classifier that generalizes to new subjects is quite challenging and in 
most cases the classifier performance can approach the subject dependent classification perfor- 
mance. It should be observed that these studies analyze EEG signal for BCI whereas in our 
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study our goal is to minimize the variability problem in fNIRS based BCI systems. 

To our knowledge there exists no study focusing on the inter-subject and inter-session related 
variability for fNIRS based BCI. Most fNIRS based BCI studies follow the research paradigm 
of developing subject /session dependent classifiers for a particular BCI task. In this study, we 
propose a classification method that is able to distinguish patterns regardless of subject-invariant 
features and subject-specific features available in the analyzed neural data. The method allows 
one to approximate a unique function in order to decode multiple recordings of fNIRS signals 
and estimate the subject-invariant feature patterns, hence addressing the issue of inter-subject 
variability, while reducing the need for the development of subject/session dependent classifier 
modelling for the subject-specific parameters. The type of functions we analyze are based on 
statistical learning theory as proposed by Vapnik [19] for pattern recognition purposes. In 
particular, we present an SVM learning framework and assess the quality of an estimate by a 
risk function R(a) and regularization framework defined in [20]. 

2. Materials and Methods 

2.1. Experimental paradigm 

We used a multi-channel optical brain-function imaging system fNIR-300 for data acquisition 
(BIOPAC Systems Inc., USA). The system has a flexible sensor consisting of 4 light sources 
having 3 built in laser diodes emitting different wavelengths of 730 nm, 805 nm, and 850 
nm and 10 detectors designed to image cortical areas underneath the forehead. With a fixed 
source-detector separation of 2.5 cm, this configuration generated a total of 16 signal channel 
measurements with a sampling rate of 2 Hz. We measured concentration levels of Oxy-Hb, 
deOxy-Hb, and total-Hb from seven healthy right-handed subjects (age range yrs old), of whom 
four subjects had no previous experience in mental imagery tasks. In this study, our focus was 
on the analysis of Oxy-Hb signals, because the changes in the concentration levels of oxy-Hb are 
highly correlated with the regional cerebral blood flow (fCBF) [21], and an increase in rCBF 
reflects an increase in neural activity [22]. The optical fiber probes were placed on the pre- 
frontal cortex of the participants following the 10-20 international electrode placement system 
[23]. Figure 2 (a) shows the locations of the transmitter laser diodes and receiver optodes which 
correspond to the left frontal region (Fpl-Ch8, AF3-Ch7, AF7-Ch6, Fl-Ch5, F3-Ch3, F5-Chl, 
F7-Ch4, F9-Ch2) and to the right frontal region (Fp2-Chl0, AF4-Ch9, AF8-Chl2, F2-Chll, 
F4-Chl3, F6-Chl5, F8-Chl4, F10-Chl6) of the brain cortex. The red-colored numbers represent 
the transmitters, and the blue-colored numbers represent the receivers; the yellow circles show 
the locations of the 16 channels. 

We divided the data recording into eight sessions with one session per day. Four of the eight 
sessions were set for performing the task of [rest — > right imagery task] and the remaining four 
sessions were organized for [rest — > left imagery task]. Thus each session was divided into two 
blocks. The duration of each block in a session was set to 3 minutes, and, overall lasted 12 
minutes for one subject. During the rest block, participants were asked to sit and relax and not 
perform any task; the signals obtained during this relaxation period were used as a reference 
in the classification task. In task blocks, participants were instructed to perform an imagery 
movement task of their right forearm in the right and left directions. This mental task was 
performed simultaneously with a given visual observation task. The visual observation task 
consisted of an action movie clip showing the movement of the participants forearms in the 
intended direction. Each movement in the clip was repeated every two seconds, which made up 
a total of 90 consecutive repetitions of the forearm movement during the task block. Figure 2 
(b) provides more details on the data acquisition paradigm. 

The present experimental paradigm appears rather different from well known motor imagery 
or cue-based BCI paradigms. Here, we propose another approach with the hypothesis that the 
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Figure 1: Experimental set-up for fNIRS data acquisition: (a) locations of the fNIRS optodes and channels in 
pre- frontal brain cortex according to 10-20 electrode placement system (inter-optode distance = 25 mm); (b) 
experimental protocol for data acquisition consisting of rest-task-rest segments; task segments last 180 s equally; 
each rest segment lasts 180 s. Subjects were instructed to perform an action observation task by simultaneously 
focusing on a video-clip that displays the movement of subjects forearm in right/left directions. 



putative human mirror neuron systems, which predict and interpret ones own actions and the 
actions of others, could be exploited for fNIRS-BCI systems, as was done in [24]. The motivation 
for this approach is based on several of the earlier mentioned studies [25-26], which successfully 
showed the feasibility of the BCI task. We verify our experiment and hypothesis using our 
test bed system, which includes the fNIRS system and the haptic devices. The communication 
among subject-fNIRS-Robot is obtained through TCP/IP protocol. Our target control system 
is the PHANTOM 1.0 premium with haptic device (19.5 cm 27 cm 37.5 cm workspace size, 
0.03 mm 3-D positional resolution increment, and 2 active degrees of freedom (DOF)). In the 
following sections we describe the proposed algorithms suitable for accurate decoding of task 
relevant fNIRS signals. 



3. Problem Formulation 

We denote x £ j^{mxn} ag a neura j signal set associated with {p«}J =1 participants performing 
{tj}j =1 cognitive tasks from {sfc}f =1 sessions; We want to estimate a unique function f(x,a*) 
from the set of functions {f(x, a), a € A} that makes similar predictions from {x} across varying 
{p, s}. Let us consider the quality of an estimate of f(x,a) as R^ pi Sl tl y(x,a) obtained for 
subject {pi} and there is another estimate R{ P2 ,s 1 ,t 1 }( x i a ) f° r subject {^2} both participating 
in a session {si} and performing the same cognitive task {tl}. Then our ultimate task is to 
minimize the following variance of both estimates 

min(V) : V = |-R{ p i, s i,ti} - R{p2,sl,tl}\, (1) 

by finding the optimal parameter {a* £ A} of a unique function. As a result the opti- 
mized function will provide statistical stability and robustness regardless of the inter-subject 
variability and inter-session variability that arise in signals acquired from different subjects 
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and sessions. We define our problem formulation in a supervised learning framework, that is, 
given {(xi,yi)f =1 }, Xi G X, an observation y G {±l} n , and a prediction f(xi,a*) denoted by a 
triplet (x, y, f(x, a)) G X x y x y. We estimate a function / : X — > {±1} which tells us whether 
a new signal observation x G X belongs to class +1 or —1. We defined the y G {±1}^ values 
for mental tasks [fa, +1), (6, —1)] as patterns related to right movement (y = +1) and baseline 
(y = — 1). Similarly, we defined [fa, +1), {b, —1)] as patterns related to left movement and base- 
line (y = —1). We perform pairwise classification in which for each possible pair of classes we 
train a separate classifier function. For M-classes then we train (M — l)M/2 classifier functions; 
thus each function is modeled to discriminate between imagery right movement and baseline 
task and for discrimination between imagery left movement and baseline. As mentioned earlier, 
our goal is to approximate a function f(xi,a), Vx G X in order to assess the quality of an 
estimate or a risk function R{ p s t y(x, a)- defined by the following structural risk minimization, 



R(x, a) = — 



£(y* ~ f(xi,a*)) 



+ AQK). (2) 



The first term in Eq.2 computes is the empirical loss function R emp [f]- which incurs a 
loss R emp [f] = for correctly classified x- otherwise R em p[f] G (0,1]- The second term is 
a regularization term defined by the Vapnik-Chervonenkis dimension which allows to control 
the capacity of the decision function, (i.e. generalization ability, convergence or confidence 
interval). By minimizing both terms, one achieves a global minimization of the risk function 
R(x,a). Thus our goal is to find some rules which allow us to minimize both terms. The 
structural risk minimization approach invented by Vladimir Vapnik leads to the development 
of support vector machines and other kernel based learning methods. 



4. Support Vector Learning 

Developed from statistical learning theory [27], the standard SVM classifier has achieved 
superior performances in almost all BCI classification tasks [28] . Because of its perfect general- 
ization ability and its well established theoretical foundations (using regularization theory and 
convex optimization) it has become the first choice in many pattern recognition tasks. The SVM 
theory has resulted in a number of extensions, such as the elegant zv-SVM [29] or leave-one-out 
machines [30] and further extensions to the case of multiple kernel learning [31]. The training 
of standard SVM requires the solution of a quadratic programming (QP) optimization which 
sometimes may not be easy to implement, in particular for large-scale problems. Therefore, we 
model SVM classifiers that are based on linear programming optimization methods (denoted by 
LP-SVM). These type of LP-SVMs are found to be sparse (kernel) classifiers efficient for very 
large scale data samples. The LP-SVM was first developed by Vapnik et al. as an extension of 
QP-SVMs; the details are provided in [19]. 

The solution of LP-SVM is based on the following linear programming optimization, 
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aeR n ,b&R I C f-f n f- \ 
i=i t=\ 

subject to yA ^a j y j K(x i ,z j ) + b\ > 1 - & 



&>0,Vi,j = l,2,...,n. (3) 

Here, a = (a±, ...,a n ), £j are slack variables and a user defined regularization parameter 
C = C(n) > 0. The corresponding decision function (a classifier) is given by, 
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Figure 2: Geometric interpretation of SVM and its decision boundary for a given toy dataset in 2D. The support 
vector set [cfi(zi), (^(22)] and forms the decision boundary with a margin M, containing all the information 

necessary to predict the class of any random test patterns. 



f(x,a,b) = sign 



(4) 



Both in Eq.2 and in Eq.3, the K(x,z) is called a reproducing kernel function endowed with 
a dot product k : X x X — > R over the elements of input patterns; this gives rise to a positive 
definite Gram matrix Kij := (xi,Zj), K E j{ ny - n [20]. In fact, this matrix contains all the 
information in order to perform data analysis and modeling of a kernel based algorithm. The 
elements of K{j corresponds to dot products of analyzed patterns in a possibly high-dimensional 
dot product space % according to the reproducing kernel property, 



K (x , z) = (<j)(x) , 4>(z)} , for any x, z G X. (5) 

This property states that 3 a feature map x — > (j)(x) and $ : X — >• H such that the evaluations 
of the kernel at points x and z are equivalent to taking the dot product between the images 
4>{x) and <p(z) of maps in some feature space H. The feature space H is known as a reproducing 
kernel Hilbert space (RKHS). An important property of the RKHS is that any Cauchy sequence 
in this space converges to a limit, in other words, it is possible to approximate a solution (e.g. a 
point with maximum similarity) as accurately as needed. Another important property of RKHS 
is that the map $ : X — > % reduces the non-linearly separable data of the original input space 
X to linear data in T~i. In view of the fact that it is only through the kernel matrix that the 
SVM receives information about the feature space and input data, we will take a close look 
at the structure of the feature space in which our data exists. The next section presents the 
basic geometry of the data in RKHS, which will lead to useful results and that can be solely 
computed in terms of the kernel. Overall, we aim to improve the BCI decoding performance by 
minimizing the inter-subject and inter-session variability. 



4-1. Kernel Feature Space 

Let us consider a finite subset S = {x\, x n } of neural data set X, a kernel function k(x, z) 
and a feature map </> into the feature space % satisfying Eq. 4. Let <&{S) = {4>{x\), (p(x n )} 
be the image of {S} under the map 0, hence (j)(S) C Ti. We investigate for the information 
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that can be obtained about $(S) in terms of its characteristics related to inter-subject/session 
variability all in RKHS. 

In RKHS, the set 3>(5) € R d forms the smallest convex polyhedron co<J>(S) that contains all 
neural data elements, 

n 

co$(S) := {<f>(x) | </>(x) = <Xi<t>(xi)}, (6) 

i=i 

where n £ N, a« > and J2?=i a i = 1- Recall that SVMs find the maximal margin 
separation in RKHS between two classes of a data set. As an example, suppose the given 
two training sets [3>(SL.),y = —1] and [<J>(S+),y = +1] together with targets and they are 
linearly separable. Thus the corresponding convex hull of both classes do not intersect i.e. 
co(3>(S_)) n co(3>(S + )) = 0, (cf. Fig. 2). Then the SVM decision boundary is formed by 
maximizing the margin between few nearest neighbor training points <p(x\) € co(3>(»S_)) and 
(j)(zi), <f>{z2) € co{<&(S + )) touching hyperplanes [H-,H + ]. Note that each 4>{x) € &(S-) for 
which <j){x) co(<J>(S_\{(/>(x)}) is called a vertex of the co($(5_) or an extremal point of the 
set. Those points are particularly interesting to us, because they strongly influence the 
forming of the decision boundaries between the two classes of the data set. Moreover, we should 
mention that inter-subject variability is closely related to different forms of decision boundaries 
and the locations of the convex hull vertices. In SVMs, the nearest neighbor vertexes between 
two convex hulls are called the support vectors, and are found as the solution of the optimization 
problem given in Eq.2. Formally stated, the Karush-Kuhn- Trucker (KKT) complementary 
slackness conditions states [32] that if a,b,£ are optimal solutions to the dual in Eq.3, then 

Vi ( S a jVj K ( x i> z j) + b ) ~ 1 + & = (7) 

M=i ' 

Mi = 1 , . . . , n 

Therefore, for any training point Xj, either a = or Ui(J27=i a jUjK( x ii x j) + b) — 1 + £i = 0. 
Because the points for which a = do not appear in the solution of the classifier in Eq.3, 
the classifier is based only on those points x« (support vectors) for which aj > that define 
the decision boundary. Therefore, these are the only points required for classification of any 
test point z. Intuitively, the SVM defines the decision boundary of convex hulls based on a 
weighted sum over the training points as in Eq.6. Those for which the weight is zero are far 
enough inside the decision boundary defined by other points to incur no loss. In other words, 
the points that are far away from the convex hulls of the data are irrelevant and will not improve 
the classifier generalization error. Notice that, in order to classify a new test point z, we only 
need to calculate the value of K(xi,z) when a > 0. Therefore, once an optimal a is learned, 
all of the training points X{ where a = can be discarded. It turns out that support vectors 
are also subject-dependent, i.e. because of variability in the data different subject-dependent 
feature space regions are formed. Then, a correct classification of any new test point depends 
on support vectors. Now, we would like to turn our attention to the role of the kernel function 
K(xi,z) in constructing the decision boundary. When, K is a simple linear kernel function 
k(xi,z) = (xi,z) we end up finding only the optimal a because the kernel does not have an 
additional hyper-parameter to be optimized. On the other hand, if the kernel is defined as 
a polynomial kernel k(xi,z) = ((xi,z) + c) p one is left to finding the optimal degree of p for 
mapping in together with a. Our analysis is based on one of the most popular kernel functions 
used in practice the Gaussian radial basis function kernel. It is defined as 
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Figure 3: Illustration of different decision boundaries in kernel feature space formed with varying Gaussian kernel 
parameters a and support vector sets; A small value of a forms a complex decision boundary with many SVs on 
the vertices of the convex hull of the data as shown in (a) a = 0.4 and (b) a = 0.6. Further, a large value of a 
will give a smoother decision surface and more regular decision boundary; thus only a few support vectors will 
have strong influences over a larger area, as shown in (c) - a = 1.2 and (d) - a — 2.4. 



This kernel measures the Euclidean distance between x% and Zj data points. In SVM, it is 
equal to measuring the distance between support vectors X{ and a random testing point Zj. In 
the feature space constructed by a Gaussian kernel, the support vectors will be in the center of 
the bell-shaped bump (e.g. Gaussian) and the a-value will determine the area of influence of a 
particular support vector over the feature space. Consequently, depending on the value of a a 
decision surface in T~L will change its decision boundary. For example, if we use a large value of a 
it will result in a smoother decision surface with a more regular decision boundary. As a result, 
a Gaussian kernel with large a will allow a support vector to have a strong influence over a 
larger area. In addition, a larger a value also increases the a value (the Lagrange multiplier) for 
the classifier in Eq.3. Also, when one support vector influences a larger area, all other support 
vectors in the area will increase by an a. value to counter this influence. Hence all o-values will 
reach a balance at a larger magnitude. It should be observed that in Fig. 2 (c, d) a large a- value 
will reduce the number of support vectors, since each support vector can cover a large space 
when only fewer are needed to define a boundary. It becomes obvious that finding both a and 
its optimal hyper-parameter a is of great importance in order to construct a decision surface in 
RKHS. Then the classification accuracy of an SVM for any random point is known to depend on 
some similarity measure (e.g. distance) to the support vectors and the decision boundary that is 
found with optimal a, a. Notice that up to now we have discussed about a single feature space 
and its decision boundary with optimal parameters. To show the BCI inter-subject variability 
in RKHS, let us now assume that Fig. 2 plots the subject-dependent feature spaces from four 
subjects when performing the same task. Then when a test point is provided to a classifier, it's 
prediction is based on the type of feature space it is analyzing the data with. Then, because of 
the strong variation in the decision boundaries a test pattern may not be classified as correctly 
as we had expected. We will present some details on the subject or session dependent feature 
spaces in the subsequent section. 
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4-2. Inter-subject variability in RKHS 

Let us now investigate the problem of possible inter-subject variability in RKHS. We de- 
note {pi,P2, ■■■,Pn} as a set of subject specific data with a set of associated feature maps 
($i(pi), $2(^2); ■■■,&n(Pn)) and optimal {o~f =1 ,a} values, where each space forms varying con- 
vex hull decision boundaries. Suppose that <&(pi) associated with {si,t\} is induced using a 
Gaussian kernel with optimal cr-value. Now recall the case, of inter-subject variability, when 
a new pattern of $(^2) associated with the same {si,ti} was provided; in such a situation 
the classifier fails to estimate the correct class for the pattern. That is because the decision 
boundary formed by cr-width for support vectors a was not appropriate. On the other hand, a 
classifier would be able to find its class correctly if we had a space $(^2) formed with optimal 
{a, a}. More formally, we introduce the following remark. 

Remark 1.1: Given an optimized function f(x,a) and a kernel k(x, •), where all x £ p\ and 
s\ with unique a; the estimation accuracy of an f(x, a) for randomly selected test point z from 
{pi, ...,p n } and {s } can be described in three different cases: 

1. Subject-dependent data: Test point : z € 4>{Pi) from {s±} were predicted accurately. 

2. Inter-session variability: Test points: z <G <p(pi) from {^2,53} were not accurately pre- 
dicted. 

3. Inter-subject variability: Test points: z € 4>{p2) , <t>{Pz) from {s\} were not accurately 
predicted. 

In the first case, one achieves an acceptable accuracy of a classifier during the test phase; 
because the testing data points belongs to the same subject {pi}, the session {s±} and the 
cognitive task {t±}. In the second case, the classifier fails to predict the class of test points even 
through they are provided from the same subject {p±} but with inter-session variability, {s2, S3}. 
The neural signal measurements obtained during the different sessions vary dramatically even 
when the same subject is used. In the last case, the classifier fails to predict the correct output 
because of the inter-subject variability problem, although the session and task were the same as 
in the first case. The latter two cases are the major problem in BCI research (i.e. the existence 
of a large degree of variability in the neural data). Every subject has different brain signals; 
also, the performance of a subject varies widely between sessions, within session, between on-line 
and off-line settings, and even from epoch to epoch. We attempt to solve these two problems 
by finding global unique features of neural signals regardless of inter-session and inter-subject 
variability. Let us consider that in case one, we succeeded in achieving accurate prediction of 
any test point as long as it belonged to the same subject and the same session. Our idea is 
based on building subject-dependent or session-dependent kernel feature space with a set of 
optimized parameters. At first glance, we might think of optimizing several functions for each 
feature space and end up with a subject/session dependent classifier. However, we want to 
avoid designing a class of functions according to the subject or the session; rather, we want to 
find a unique function that will find the decision boundary for all subjects or sessions. In order 
to avoid such problems, we consider constructing a feature space with several kernel functions 
and related optimal parameters. As an example, consider several feature spaces induced by a 
Gaussian kernel in which each space is optimized with good value according to the data of the 
subject. Out task is to incorporate all feature spaces into one richer feature space and learn the 
data with a unique SVM function. More formally, we define several subject dependent feature 
spaces that result in a subject dependent decision boundary in RKHS. 

4-3. Extension to multiple RKHS 

Since the kernel functions in Eq.5 induce the necessary RKHS for our analysis, we must 
therefore start this section by presenting some closure properties of kernels that allow us to 
create extended RKHSs. 
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Proposition 1. Let K\ and K 2 be arbitrary positive definite kernels over X x X where 
X € R is a non-empty set: 

The set of positive definite kernels is a closed convex cone, that is, (a) if a\, a 2 > 0, then 
ct\K\ + a 2 K 2 is a positive definite; and (b) if K(x, z) := lim n ^ OQ K(x, z) exists for all x, z, 
then K is a positive definite; 

The above statement is fairly obvious that a linear combination of positive definite kernels 
preserves the finitely positive definiteness of the "kernel" property. For instance, recall that a 
real symmetric n x n matrix K satisfying J27j c i°j^i,j — f° r an y vector c € R n is a positive 
definite matrix. Subsequently, for any vector c G R n , the kernel K is positive definite V c as 
shown in 

c'(^i + K 2 )c = c'K lC + c'K 2 c > 0. (9) 

The feature spaces associated with a set of multiple kernels can be defined by their linear 
combinations K = J2n=i^2i a i^n (where a £ R), by forming the concatenation of the corre- 
sponding feature vectors from all kernels. As an example, from the definition of the reproducing 
property in Eq.5, one may easily verify the following equality, 



Ki + K 2 = (0! (*)0i (*)) + (<k{x)<h{z)) 

= (<Kx),<Kz)) = K(x,z); (10) 

In general, to construct feature spaces consisting of different classes of Gaussian kernel 
functions with optimal kernel parameter at different points in the space, one can consider the 
following relationship, 

k k 11 _ i|2 

^K n {x h z) = ^<exp( — X ^ 2 Z ) (11) 

n=l n=l n 

Here, each K n induces subject-dependent optimized features spaces with the range of optimal 
parameters {(ajai), (at,<7i), (a^Cfc)}) overall giving an improved decision surface in RKHS. 
Using the last equation (Eq.ll), we are required to minimize the following optimization problem 
(see Appendix 1 for derivation), 



mm 

a£R n ,£,£R I C 



{-. k n -. n 

n=l t=l i=l ' 



subject to constraints, 



^n=U=l ' 

<>0, i,j = 1,2,..,/, n=l,2,..,/c 

^ >0,Vi = 1,2,..,/. (13) 

with respect to af and ^ and b and for which we denote af as the weight of the support 
vectors for the i — th training vector for the n — th kernel. The corresponding decision function 
using multiple kernels is given by, 
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r k I 



f(x,a,b) = sign 



• n=l i=l 
.,2„.2 ^ ^_ , k„,k , 



... + atyiK 2 ( Xi , z) + afyfK k ( Xi , z)) 



(14) 



where it works with multiple sets of kernels by choosing the values of af and <Ji, i = 1, .., n, 
rather than just values and a single er-value as given in Eq. 4. When using Gaussian kernel 
functions of different widths (<r- values), one then approximates the unknown function with a 
set of decision functions of the form: 



f(x,a,b) = sign 



n=l i=l 



2al 



(15) 



The function in Eq.15 allows us to choose a decision rule from a suite of different classes 
of function, and to choose the optimal kernel parameter at different points in the space, for 
instance as shown in Fig.3. 



5. Experimental Results 

We performed our experiments on the fNIRS data sets detailed in section 2. In total, we 
performed four experiments and model classifiers according to the experiment type. The first 
two experiments were done using standard LP-SVM formulations and the remaining two were 
done using multi-kernel LP-SVM formulations. In each experiment, we re-arranged the analyzed 
data set in order to suit the experiment type. 

During the modelling of classifiers in each experiment, we conducted an optimization of the 
hyper-parameters C (regularization parameter) and a (Gaussian kernel parameter) for every 
LP-SVM classifier. As mentioned earlier, their combined values influence the decision boundary 
in RKHS, and thus affect the classification performance. In order to perform this optimization, 
different methods such as grid search or gradient descent algorithm can be used. The most 
practical approach for determining the optimal parameters is usually based on evaluating a 
cross validation error rate and the performance of the classifier model in order to minimize the 
overall test error (risk function). Here, we performed an extensive search on optimal C and 
a using 10-fold cross validation. We tried various pairs of (C, a) and selected one ,for each 
subject, with the minimum cross validation error. Algorithm 1 describes the steps that were 
done in order to model a subject specific classifier. Below, we present each experiment in more 
detail and elucidate its purpose. 

The performance levels of the classifiers are produced in terms of their cross validation 
error and generalization test error. In order to comprehend the achieved results and their 
significance from each classifier, we introduce three classifier performance regions in order to 
define best/worst and acceptable classification rates. The performance regions are defined as 
follows, 

1- Qbest ■= [0, ...,0.11] if Error rate < 11% , 

2. Qaccept ■= (0.11, ...,0.20] if 11% < Error rate < 20% , 

3. Qworst ■= (0.20, 0.50] if 20% < Error rate < 50% . 

Thus the performance of a classifier will be judged according to the above criteria. Certainly 
our first goal is to achieve the best possible classification results in the Qbest region. 
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Algorithm 1: Standard LP-SVM training 10-fold CV 
Input: {xi,yi}f =1 ; Xi G X 
foreach Experiment 1 or Experiment 2 do 

• Perform 10-fold CV and a search on (C, a) 

Divide the training set (of size) m into n disjoint 
sets Si,...,S n of equal size n/m 

foreach Si do 

• Train a LP-SVM (Eq.3 and Eq.4) on S/Si 

• Test it on Si •->■ error(i) 

Output: average test error : error(i)/n 
Output: 

• Set of optimal parameters for each experiment, 
a , a = {<7i, (72, crj}' p and C 

• Minimum generalization error 



5.1. Standard LP-SVM classification 

We performed two experiments using standard SVM modelling by considering inter-subject 
variance and inter-session variance. One purpose is to find the optimal kernel parameters, 
according to experiment type, that form the optimal decision boundary. For more detail, a 
reader is advised to refer to Algorithm 1 which provides more detail on the algorithmic flow 
corresponding to Experiment 1 and Experiment 2. 

Experiment 1 (Subject dependent classification) 

For each subject {pi}J =1 and task pair [ti,b] and [t2,b], we train and test separate subject- 
dependent classifiers in total making up 14 classifiers. Seven of these are trained for classification 
of [t\,b] and the remaining seven are for [i2,&]- The total fNIRS data were re-arranged into 
seven sets (Xp^' b \ Xp£' b \.., Ap* 1 ' 6 ') and seven sets (Ap* 2 ' 6 ', xjj£' b K.., Ap* 2 ' b '). Each set includes all 
fNIRS recording from eight sessions according to the subject type. As explained earlier, in 
order to perform cross validations we divide each set into 10-disjoint subsets. The expected 
output of this experiment consists of all subject-dependent optimized classifiers with minimum 
generalization error and the optimal classifier hyper-parameters [C, a]. This experiment was 
conducted to test and assess how well a subject specific optimized classifier performs on the 
given neural dataset. To start, let us first discuss the results obtained from experiment-1, 
provided in Table 1. Here, the first column represents the subjects' data re-arranged according 
to this experiment type. The next column, defines the type of imagery tasks (recall [t\ -left 
and ti- right]). Here and later, each second row of the table, colored in gray, represents the 
results associated with the right imagery movement tasks. The last two columns of the table 
shows the average generalization error obtained from 10-fold cross validation (i.e. the average 
of 10 different CV errors; also see algorithm 1). In the next column, we present a random 
test experiment for the optimized classifier model from the CV phase with minimum CV error. 
In this case, the number of input test samples varied X € ^mxi6 depending on the fNIRS 
recordings in which m- could be in m G [3, 10]. It is suggested that one analyze the produced 
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tables in a column-wise order for easy comparison, for instance, the performance of all subject- 
dependent classifiers ranging from [si, s7]' as well as the number of support vectors obtained 
from optimal SVM model. The optimal hyper-parameters [C, a] are seen to vary according to 
the subject specific trained/optimized classifier model. It is also a common experience that CV 
errors usually tend to be less than test errors. This is because during the CV phase the random 
test patterns are provided from the analyzed data set, while, test errors are the results for any 
random test point representing of of [ti , ^2] pairs. For every subject, the CV errors lie within 
the ®best region except for subject 4. However, no types of CV error were below the acceptable 
performance region. Most importantly, one needs to achieve a minimum in generalization test 
errors when a new unlabelled pattern is given, as shown in the last column of the table. 

Here, the best possible test error is obtained from a classifier modelled for subject 3 in 
recognizing right imagery movement (this is almost 94% accuracy). We have also noticed, 
however,that the worst classification performance was obtained for subject 4 in both cognitive 
tasks. The majority of the subject-dependent classifiers performed in the acceptable region 
and some performed in best region (bold text in Table 1). This is thought to be the optimal 
case, because each classifier was trained from the data of the same subject during the various 
sessions. Even when the data came from the same subject, it was inter-session data, and thus 
there was inter-session variability. We could conceivably be satisfied with the current results, but 
one doesnt know if the subject-dependently modeled classifier would predict any test points from 
other subjects (assuming the underlying neurological process is the same when performing the 
same cognitive task). To examine this issue, we have extensively analyzed the exact performance 
of all classifiers for all subjects in total, making up 49 experiments during which test patterns 
were given from different subjects to different classifiers. 



Table 1: Experiment 1: The performance of subject-dependent LP-SVMs classifiers in terms of their cross 
validation (CV) and test errors. Here the analyzed input fNIRS datasets are from and [t2,b]. The best 

found hyperparameters C and a for each classifier is shown in columns 3 and 4. 



Data 



Task Type I SV C a CV Error Test Error 



Subject 1 left task 




Subject 2 


left task 


185 


100 


9.34 


0.10 


0.15 








Subject 3 


left rask 


156 


100 


0.84 


0.09 


0.18 




right task 


185 


100 


0.85 


0.02 


0.06 


Subject 4 


left task 


185 


80 


1.05 


0.17 


0.24 




right task 


165 


80 


1.92 


0.19 


0.26 


Subject 5 


left task 


248 


100 


0.51 


0.10 


0.21 








Subject 6 


left task 


148 


100 


0.27 


0.11 


0.31 








Subject 7 


left task 


155 


100 


1.55 


0.08 


0.29 




right task 


150 


100 


2.64 


0.06 





Table 2 shows the results when an optimized classifier was provided with varying test pat- 
terns (first column : subject data). In each column one may see the performance of each 
classifiers test error. Notice that the diagonal entries of the table represent the results obtained 
when the input test patterns were provided from the same subject (similar to our previous 
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experiment). Under our three performance regions, if we interpret the table, one can notice 
that the diagonal entries of the table never go below the acceptable region; moreover, 60% of 
the obtained results are in the best performance region. Nevertheless, the other column-wise 
entries result mostly in worrisome results, in which many error rates fall below the acceptable 
region 

The worst classification test error was 28% in [fa, subject 4]; the best result when considering 
inter-subject variability is in [f 4, subject 1]. We have bolded the text in the table (the error 
rates € Qbest) to make it easily distinguishable from the remaining worst performances. In 
brief, we can say that classifiers modelled in a subject-dependent manner do not perform well 
when tested against other subjects' data. This problem again clearly shows the presence of 
inter-subject variability in the data. 



Table 2: The performance of subject-dependently optimized classifiers when input test patterns were provided 
randomly and independent of subject type. Here, {/i}J=i represents the seven different classifiers trained for each 
subject The diagonal entries of the table (when i = j) represent subject-dependent classification (same 

as in Table f). Other entries represent the classification test errors when considering inter-subject variability. 



Data 



Subject 1 



Subject 2 



Subject 3 



Subject 4 



Subject 5 



Subject 6 



Subject 7 



fi 



h 



h 



h 



6 



f 



7 



0.2145 0.2750 0.0925 0.0320 0.1810 0.1440 0.2010 
0.0650 0.1400 0.0785 0.2100 0.2100 0.2140 0.4655 



0.2150 
0.1895 



0.0320 
0.0810 



0.1931 
0.2265 



0.1460 
0.1125 



0.1715 
0.1405 



0.2360 
0.1715 



0.1230 
0.1440 



0.2801 
0.1345 



0.2312 
0.1909 



0.1421 
0.1561 



0.1901 
0.2511 



0.2311 
0.2205 



0.2603 
0.1571 



0.1303 
0.1640 



0.1502 
0.1952 



0.2311 
0.1815 



0.309 
0.2161 



0.1101 
0.1225 



0.2512 
0.2035 



0.2612 
0.2125 



0.2830 
0.1140 



0.1112 

0.1951 



0.1320 
0.1102 



0.1201 
0.1651 



0.1461 
0.1315 



0.1575 
0.2403 



0.1209 
0.1811 



0.1350 
0.1740 



0.1511 
0.1091 



0.1521 
0.1811 



0.1921 
0.1285 



0.1110 
0.1025 



0.1215 
0.1210 



0.1361 
0.1105 



0.1221 
0.1240 



0.1450 
0.1204 



0.0620 

0.1320 



0.1912 
0.1625 



0.1460 
0.1121 



0.1715 
0.2105 



0.2312 
0.1735 



0.1230 
0.1340 



In our next experiment, we analyze the inter-session variability problem that arises in a 
manner similar to that of inter-subject variability. 

Experiment 2 (session dependent classification) 

In this experiment, we model session-dependent classifiers in the same manner as in Ex- 
periment 1. The difference is that here the input data were re-arranged so that each set rep- 
resents each session from all eight sessions separately. Thus, in total, we have eight datasets 



M , Xll 1 X£ ' 0| ) for [t ! , b] and eight datasets ( ' 0| , X^ ' OJ . . . , X^ ' OJ ) for [t 2 , b] . Here, 
each set consists of NIRS data recorded from all seven subjects according to the session number. 
Thus, we train and test sixteen session-dependent classifiers. As the output, we obtain opti- 
mized classifiers with average generalization error for each session and optimal hyper-parameters 
[C, a]. This experiment was conducted to test and estimate how well a session specific optimized 
classifier performs on the dataset. 

Table 3 presents the overall classification performance in terms of CV error, and test errors 
together with the found optimal classifier hyper-parameters and support vectors. The discussion 



,[ti,b] 



[*2,&] y [t%,b] 



4*2 M n 
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here will be similar to our previous discussions of Table 1, with the difference that we are here 
considering the session-dependently trained classifier outputs. In a nutshell, the performance of 
session-dependent classifiers were not worse than that of subject-dependent classifiers (shown 
in Table 1). In this experiment, we were able to obtain minimum possible cross-validation 
errors, which is a large improvement compared to results shown in Table 1. In addition, the 
test results were in the acceptable region Qaccept, with one exception which belonged to the 
best classification region. Generally, when we performed 10 random trials of the test patterns 
with the classifier, it was noticed that each time the error rate did not fall to the worst clas- 
sification region. Moreover, we have performed additional experiments in order to analyze the 
performance of the session-dependent modeled classifiers when the test patterns were randomly 
drawn from different sessions. 

Table 3: Experiment 2: The performance of session-dependent LP-SVM classifiers in terms of their cross 
validation (CV) and test errors. Here the analyzed input fNIRS datasets are from [ti , 6] , and [t2,6]-The best 
found hyperparameters C and a for each classifier are shown in columns 3 and 4. 



Data 


Task Type 


SV 


C 


a 


CV Error 


Test Error 


Session 1 


left task 


148 


150 


3.88 


0.02 


0.28 
















Session 2 


left task 


125 


100 


4.66 


0.07 


0.13 




right task 


135 


100 


3.24 


0.18 


0.24 


Session 3 


left task 


118 


100 


3.22 


0.14 


0.18 




right task 


168 


100 


2.64 


0.00 


0.16 


Session 4 


left task 


178 


80 


1.39 


0.06 


0.17 




right task 


154 


80 


2.09 


0.10 


0.23 



Table 4 shows the experimental results obtained across session-independent random test 
patterns as inputs to classifiers modeled in a session-dependent way. The diagonal entries are 
assumed to define the best performance, because the test patterns come from the same trained 
source. The problem with this approach is that the modeled classifiers mostly fail to predict the 
correct outputs for test patterns from irrelevant sessions. Let us for example simply compare the 
off-diagonal entries to the diagonal elements of the same table. Then, we can notice how much 
the performance of the classifiers degrades, dropping out of the acceptable region. We calculated 
the average error rate of off-diagonal elements in column- wise order that did belong to the Q WO rse 
region. Such a situation can be seen in the results in Table 2 as well. We have established such 
non-robust classifier performance through several trials (most importantly empirically) . In the 
following Experiments we show how one can overcome this problem or reach the Qaccept or 
@best region by modeling a multi- kernel LP-SVM approach (explained in section 5). 

5.2. Multi-kernel LP-SVM classification 

Recall that our aim is to minimize inter-subject or -session variability by improving the 
classifier performances provided previously in Table 2 and Table 4. The following experiment 
was performed using the LP-SVM formations provided in section 5, namely, the SVM that 
learns multiple RKHS induced by the Gaussian Kernel and according to the subject type. We 
mentioned that it was possible to construct subject dependent feature spaces (RKHS) and to 
consider the convex unions of several such spaces in order to find global unique neural signal 
features in RKHS. This assumption is supported empirically in the following experiment types. 
More details on the algorithmic flow of the upcoming experiments are provided in Algorithm 2. 
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Algorithm 2: Multi-kernel LP-SVM training 
Input: {x?,yi}f =1 ; a = {01, cr 2 , a 7 } and C 
foreach Experiment 3 or Experiment 4 do 

• Perform multi RKHS mapping x $i(x) = (<f>n(x), 4>iki ( x )) € ^f' 
with the range of parameters a and C using 

K n (x,z) = exp(-||x - z\\ 2 /(2al), n = 1.,,,.5 

• Perform 10-fold cross validation 

Divide the training set (of size m) into n disjoint 
sets S\,...,S n of equal size n/m 

foreach Si do 

• Train a multi-kernel LP-SVM (Eq.14) on S/Si 

• Find optimal a* by solving Eq.12 and Eq.13 

• Test the optimized classifier 

• Estimate / : Si i-> {±1} 

_ Output: average test error: error(i)/n 
Output: 

• Optimal multi-kernel LP-SVM model 

• Minimum generalization error 
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Table 4: The performance of session-dependently optimized classifiers when input test patterns were provided 
randomly and independent of session type. Here, {/i}* =1 represents the four different classifiers trained for each 
subject {pj}1- =1 . The diagonal entries of the table (when i = j) represents session-dependent classification (same 
as in Table 3). Other entries represent the classification test errors when considering inter-session variability. 



Data 


h 


h 


h 


h 


Session 


1 


0.1711 


0.1932 


0.1232 


0.0920 






0.2501 


0.2312 


0.1319 


0.2891 


Session 


2 


0.2312 


0.1872 


0.0982 


0.2012 






0.1672 


0.1298 


0.2130 


0.3645 


Session 


3 


0.2321 


0.1892 


0.1011 


0.2509 






0.1912 


0.2310 


0.1982 


0.2490 


Session 


4 


0.1710 


0.2117 


0.1821 


0.1801 






0.2019 


0.1789 


0.0951 


0.2100 



Experiment 3 (subject independent classification) 

In this experiment, we define subject-dependent multiple-kernel feature space construction 
for the fNIRS data set explained in Experiment 1. We define seven Gaussian kernels with a the 
range of a parameters for (Xp^' b \ Xp£' b K.., A^* 1 '^) and seven kernels for (Xp?' b \ Xp£' b K.., Xp 7 2 ' b ^). 
We further perform parallel learning of two multi-kernel LP-SVM classifiers (as described in 
section 5, Eq.16), one for classification of [ti,b] and the other for [i2>&]- Notice that in Experi- 
ment 1 and Experiment 2 we optimized multiple subject /session dependent classifiers. As shown 
in Algorithm 2, after defining the input dataset with the range of optimal hyper-parameters 
[C, a], we perform kernel mapping. In other words, we define the subject dependent feature 
spaces each optimized with support vectors and classifier parameters such as a, a etc. Further, 
we perform 10-fold CV in order to find the optimal decision surface made of with the support 
vectors and the related parameters. In this experiment, the cross validation error rate also 
plays an important role, because by combining multiple feature spaces, we are re-defining the 
decision boundary when the input patterns consist of random subject-independent data sets. 
Let us now consider the outcomes from this experiment, as shown in Figure 4. The (A)-plane of 
the figure demonstrates the generalization error at each iteration through the 10-fold CV. For 
a left imagery task [t\, b] the minimum CV error rate was obtained from 7-10 fold in the range 
of [0.11, 0.11, 0.10, 0.09]. Similarly, in a right imagery task \p2,b\ the classifier performance 
was around [0.04, 0.11, 0.06, 0.07]; which indicated better classifiable performance. In both 
cases, after the 3-fold the CV error rate drops into to the best performance region. Further, 
the (B)-plane in Figure 4 shows the test results for randomly drawn patterns from all seven 
subjects as indicated in the x-axis. Here, the obtained test errors are taken from 5-trials of 
random input from each subject over all subjects. For almost all subjects the test error rate 
appears in the best performance region @best- This is except for one case for the test patterns 
from subject three (right imagery task), for which the error rate reaches 0.16 ± 0.3; However, 
this is still in the Q acC ept region. To see the significance of these results, we compare them to the 
data in Table 2, the results of which were obtained from subject-dependent classifiers. Overall, 
performances of the multi-kernel SVMs for any randomly taken test patterns are superior to the 
off diagonal performances, shown in Table 2. Another important advantage of the latter method 
is that, unlike training multiple-subject dependent classifiers, we train only a single classifier. 
The classifier is empirically shown to perform better even when the input test patterns are 
independently taken across varying subjects. Moreover, compare the total number of support 
vectors (making the decision boundary for multi-kernel SVM) and the number of support vec- 
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1 2 3 4 5 6 7 8 
(A) - CV trials (10-fold) 



9 10 



0.35 
0.30 
0.25 



£ 0.20f 



0.15- 
0.10 
0.05- 




-[H,b] 
-[t2,b] 




1 2 3 4 5 6 7 
(B) - No. # Subjects 



Figure 4: Subject-independent classification results using multi-kernel LP-SVMS. (a) - tenfold cross- 
validation error rate for both [ti,b] and [tijia] tasks. Here the classifier is trained by using the dataset from 
all subjects ; (b) Performance of the optimized classifier from (a) in terms of its test error rate with input test 
patterns drawn randomly from subject-independent dataset. 




1 2 3 4 5 6 7 8 
(A)-CV trials (10-fold) 



9 10 



0.30 
0.25 

£ 
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o 

LU 0.15- 
0.10 
0.05- 

L 



-[t1,b] 

-[t2.b] 




2 3 
(B) - No. # Sessions 



Figure 5: Session-independent classification results using multi-kernel LP-SVMS. (a) - tenfold cross- 
validation error rate for both [t\,b] and [ti,fe] tasks. Here the classifier is trained by using the dataset from 
all sessions ; (b) Performance of the optimized classifier from (a) in terms of its test error rate with input test 
patterns drawn randomly from session- independent dataset. 



tors with subject-dependent classifiers. The obvious thing is that in multi-kernel representation 
the number of support vectors increases. This is due to the fact that while learning from a 
dataset with variance the optimal subject-dependent decision boundaries merge or adapt to the 
whole training set. In the subsequent experiment, we will analyze the data with inter-session 
variability using the same approach above. 

Experiment 4 (session independent classification) 

Our final experiment was focused on minimizing the inter-session variability. It is performed 
in a similar paradigm as that of Experiment 3 except for there being different input data set 
arrangements and optimal hyper-parameters. We define session-dependent multiple-kernel fea- 
ture space for data (AfJ' 1 '^, xjfe ^J* 1 ' ) and (<Yj* 2 , A*]* AfJ* 2 '^) and train two parallel 
multi-kernel LP-SVMs. Clearly, we avoid training and testing multiple session-dependent clas- 
sifiers by instead using multiple kernel induced feature space and two classifiers. The empirical 
results for this experiment are provided in Fig. 5, in which the (A) - plane demonstrates the 
CV error rate across varying folds and the (B)-plane shows the results associated for a ran- 
dom test experiment. First, we should notice that this experiment outperformed the previous 
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experiment-3 during both the CV and the Test phase. Perhaps, this result is due to the lower 
number of sessions. In Fig.5(B), the test error rates are significantly improved compared to the 
data shown in Table 4. Most of the obtained results in our last experiment belong to the Qbest 
region, while the results fluctuated in experiment 2. The best performance of a classifier was 
set for any test points from subject 1 and subject 2 considering the imagery right movement 
task. Through this experiment, we have shown that it is possible to overcome the issue of BCI 
inter-session variability and to avoid modeling of the session-dependent classifiers. The resulting 
function will find patterns related to varying subjects/sessions in an adaptive manner. 

6. Conclusions 

We have presented a methodology to model a multi-kernel SVM classifier functions; This 
methodology allows us to refrain from training subject /session data dependent multiple clas- 
sifiers for specific BCI task (in standard BCI, several subject data dependent classifiers are 
trained). Instead, we incorporated all subject /session specific feature spaces into unique and 
much richer feature spaces, allowing us to form a more robust set of decision functions to be 
employed in data classification. Empirical evidence shows that our approach achieved more ro- 
bust predictions in comparison to those of standard data-dependent classification. The average 
classification test errors rates were in a 5% < ErrorRate < 10% region regardless of strong 
variations in the neural data. The effectiveness of our approach was particularly noticeable 
when the testing data were provided in our random BCI test experiments. These initial exper- 
iments have proven the validity of the new technique for a small number of subjects. However, 
the computational cost of the algorithm is proportional to the number of subjects' data. We 
used a linear-optimization approach in order to minimize the computational cost of the overall 
algorithm. Moreover, for most off-line modeling of BCI data, the computational cost is not a 
major problem. During the online classification experiments, the decision is based on the num- 
ber of the support vectors (not the whole training set), therefore this approach is computational 
feasible for real-time BCI. Our further experimental work should aim to support these initial 
results. 



7. APPENDIX 1 

In Fig. 2, we showed a feature space in R 2 in which our task was to maximize the margin 
M = jj^jj by minimizing the \\w\\ -weight vector. This problem is presented as a constrained 
optimization problem. When learning multiple feature spaces our optimization task is to con- 
sider multiple margins associated with each space. For instance, the primal SVM optimization 
form using multiple RKHS is given as a quadratic programming problem, 



mm 

w£R n ,beR 



L n=l Z i=l J 

subject to y™( E ^^(xi)) + b) > 1 - £ 



&>0,Vi = l,2,...,i; 

n= 1,2, ...,k; l,n € N (16) 

where w is the weight parameter to be minimized in order to maximize the margin (see 
Fig. 2), C and £ are the before regularization trade-off and the slack variable. The corresponding 
decision function for the multiple kernel case is defined as, 
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f(x, w, b) = sign 



n=l i=l 



(17) 



In this formulation, we are considering a set of multiple mapping functions $ n , n = 1, 2, fc 
to map data patterns into RKHS, for instance 



x ^> $ n (x) = ($!(x), $ 2 (x), $ fc (x)) € 7& 

where $j(x) £ H^, U = H® 1 x ■ ■ ■ x H* k . (18) 

We solve the primal quadratic optimization problem by finding the saddle point of the 
Lagrange function: 



1 k C ' 

L(w,b,£,a,P) = g S IWI 2 + 2"E& ~ 

ra=l i=l 
k n n 

-EE ®n(xi)) + 6) - 1 + - E 

(19) 

where a and j3 denote Lagrange multipliers, hence a« > and ft > 0. By differentiating 
with respect to u>,6 and £j, the following equations are obtained: 



^ = o => «, = 2 E «f % n <M^); (20) 

n=l i=l 

^ = => a? yf = 0; (21) 
_ = ^af = C-/?f; (22) 

Substituting Eqs.(20) and (21) into the Lagrangian function Eq.19, we obtain the dual form 
of the optimization problem, which is purely represented by dual variables, 



, k n k n n 

S EE«""EE a?a]y?yJ(<S> n ( Xi ), *„(*,-)) \ 

n=l i=l n=l i=l 

subject to a^yf = 0, 

< a? < C,V» = 1,2,...,/, (23) 

Using the above formulation, we are still dealing with quadratic programming optimization. 
In addition, we loose the sense of a multiple decision boundary defined by the multiple kernels. 
Therefore, we derive the equivalent minimization problem from Eq.16 by representing the min- 
imization problem by its dual variables a, (5 and its expansion vectors from Eq. 20-22 and using 
the Eq.5. Then, we obtain the following minimization problem using a set of multiple kernel 
combinations, 
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mm 

aeR n ,£eR 

....... %J 

subject to constraints 

k n 



-l ■ ■ -i '* -r J 



(A/ ft \ 

£ £ a?yjK n (xi, Zj) + b)>l-^ 
n=\i=\ ' 



a™ > 0, z = l,2,...,n, n = l,2,...,fe 

&>0,Vi = l,2,...,n. (24) 



This approach automatically leads to a decision function of the following form, 



f(x,a,b) = sign 



(25) 



■ k I -, 

The solution given in Eq.24 is essentially a consequence of the form of the regularizer shown 
in Eq.2 (second term) and the function in the first term of Eq.2 is given by Eq.25. 

The idea of linear programming SVMs is to use the kernel expansion as an ansatz for the 
solution, but to use a different regularizer, namely the l\ norm of the coefficient vector. The 
main motivation for this is that this regularization is known to induce sparse expansions. 

This amounts to the following linear objective function 

o EE* + ;;?>} (*> 

n=l i=l i=l ' 

k n 



ubject to yi i E E 0!iyjK n (xi, zj) + bj > 1 - & 

^n=li=l ' 



a™ > 0, i = l,2,...,n, n = l,2,...,k 

& >0,Vi = l,2,...,n. (27) 

with respect to a" n = 1,2, and £ and 6. Initially, Vapnik proposed the Linear pro- 
gramming approach for a single kernel SVM approach in [19]. He proved that for the QP-SVM, 
the expected ratio of support vectors over the size of the training data can serve as an approx- 
imate upper bound of generalization ability of SVM. His assumption is based on the idea that 
in LP-formulations, the classifier (Eq.25) has the same form of kernel expansion for the support 
vectors as do the QP-formulations; He used a linear objective function to directly minimize 
the number of support vectors of the separating hyperplane. Then still in linear programming 
formulations the separating hyperplane remains exactly in the same form of kernel expansion 
and the error constrains also remain same. Finally, our goal in using LP-SVM for multi-kernel 
induced feature spaces is then to minimize the number of support vectors in order to obtain a 
sparse classifier to analyze data in multiple RKHS. That's to say that for each feature space we 
select only those important support vectors x for a > (see Eq.7) that are crucial in making the 
decision boundary. These vectors are those points explained in Fig.2 as vertices of the convex 
hull of the data and those points in Fig.3 corresponding to different Gaussian bumps. For other 
types of sparse SVM classifiers refer to [33-35]. 
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